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Clumping in the inner winds of hot, massive stars from 
hydrodynamical line-driven instability simulations 
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ABSTRACT 

We investigate the effects of stellar limb-darkening and photospheric perturbations 
for the onset of wind structure arising from the strong, intrinsic line-deshadowing 
instability (LDI) of a line-driven stellar wind. A linear perturbation analysis shows 
that including limb-darkening reduces the stabilizing effect of the diffuse radiation, 
leading to a net instability growth rate even at the wind base. Numerical radiation- 
hydrodynamics simulations of the non-linear evolution of this instability then show 
that, in comparison with previous models assuming a uniformly bright star without 
base perturbations, wind structure now develops much closer (r < l.li?^,) to the 
photosphere. This is in much better agreement with observations of 0-type stars, 
which typically indicate the presence of strong clumping quite near the wind base. 

Key words: stars: early-type - stars: mass-loss - stars: winds, outflows - hydrody- 
namics - instabilities 
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1 INTRODUCTION 

Hot, massive stars possess strong winds driven by line scat- 
tering of the star's intense ultra-violet (UV) radiation field. 
The first quantitative description of such line driving was 



given in the seminal paper by Castor et al. ( 1975 CAK), who 
assumed a smooth, steady-state outflow. But even though 
extensions of this theory (e.g., Pauldrach et al.|1986 Friend 
|fc Abbott 1986 1 have had considerable success in explain- 
ing many global properties of OB-star winds, like the pre- 
dicted mass-loss dependence on metallicity and the relation 
between the wind-momentum and the star's luminosity, it 
is nowadays clear that these winds are in fact both highly 
variable and structured on a broad range of temporal and 
spatial scales (see Puis et al.|2008 Sundqvist et al.|2011 for 
comprehensive reviews). 

Linear stability analyses have shown ( [MacGregor et ah] 
[19791 lOwocki fc Rybicki|[T984l [TMsl the last two hereafter 
ORI, ORII) that the line driving of these winds is subject 
to a very strong intrinsic instability, operating on small spa- 
tial scales. And subsequent numerical modeling of the non- 
linear evolution of this Ime-deshadowing instability (LDI) 
has confirmed that the time-dependent wind indeed devel- 
ops a highly inhomogcneous, 'clumped', structure (Owocki 
et al.|1988][Feldincicr 1995( [Dessart fc Owocki|2003 ~ 



Such structured LDI models provide a natural explana- 
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tion for a number of observed phenomena in OB-stars, such 
Eis the soft X-ray emission and broad X-ray lines observed by 
orbiting telescopes like Chandra and xmm-newton |Feld- 

J' 



meier et al .'1997!'B erghoefer et al.|T997||Giidel fc Naze'200'9[ 
Cohen et al.,2010( ) , the extended regions of zero residual flux 
typically seen in saturated UV resonance lines ( Lucy|[T983 



Puis et aL]|1993 Sundqvist et al.|[2010 l, and the migrating 



spectral sub-peaks superimposed on broad optical recombi- 
nation lines (Evcrsberg et al. 1998 Dessart fc Owocki|2005b 



Lepine 



Mofratj^2008 1 



However, a multitude of independent observational 
studies also suggest the presence of clumps in inner wind 



regions close to the photosphere ( 


Eversberg et al.|1998 


Puis 


et al. 


20061 ISundqvist et al.|2011| 


Najarro et al.|2011||Cohen 


et al. 


2011||Bouretet al. 


2012 


1 ; this is not reproduced by con- 



servative, self-excited LDI models, which develop structure 
only away from the photospheric wind base (at r > 1.5i?t, 
[Runacres fc Owocki|[2002| . 

This result has led to a common perception that the 
LDI is not able to produce structure in inner wind regions, 
and so that some other process may be the main agent re- 
sponsible for the overall features of wind clumping in OB- 
stars. But note that the strong damping of structure in the 
inner wind found in previous instability models was a di- 
rect consequence of a complete cancellation of the LDI by 



the counteracting line-drag effect (Lucy 19841 at the stel- 
lar surface (e.g. [Owocki fc Puls|1996| hereafter OP96). But 
as already pointed out by ORII, if one accounts for limb- 
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1.2, 1.3 

r/R. 



Figure 1. Relative damping effect of the line-drag on the ORI 
pure-absorption growth rate (eqn.jsj for uniformly bright (dashed 
lines) and limb-darkened (solid lines) stellar discs. Various ap- 
proximations for the scattering source function as indicated in 
the figure, where 'FF denotes full angle integrations. Note how 
limb-darkened models are unstable also at the stellar surface 

(r2net(i?*) > 0). 



darkening of the stellar surface radiation, this cancellation 
should be incomplete and so lead to an unstable wind base. 
However, whereas the effect of such limb-darkening has been 
considered before in the context of steady-state, line-driven 
winds (Cranmer & Owocki 1995 Cure et al. 2012[ ), it has 



never been explored in time-dependent simulations. 

This paper follows up on this early conjecture regarding 
the effect of limb-darkening on the LDI. In !j2]we perform an 
analytic linear stability analysis including simple Eddington 
limb-darkening of the photosphere. i|3]then describes our nu- 
merical, self-consistent, radiation- hydrodynamical modeling 
technique for studying the non-linear evolution of the com- 
petition between the LDI and the line-drag. !|4] examines to 
what extent the new models including limb-darkening, as 
well as a simple photospheric sound wave perturbation, can 
produce significant wind structure also in the inner wind. Fi- 
nally, !j5] discusses these results, gives our conclusions, and 
outlines future work. 



2 ANALYTIC PERTURBATION ANALYSIS 
INCLUDING LIMB-DARKENING 

Let us first examine the effects of limb- darkening in the lin- 
ear regime when perturbing the radiation force exerted by a 
single line. The instability arises from perturbations in the 
direct component of the line-force, which is proportional to 
the stellar core intensity /^(/i, r) = I^D{fi,r). Here sets 
the scale for the intensity and the flux-normalized disc func- 
tion -D(/x, r) accounts for the variation of intensity along lo- 
cal direction cosine ^ at radius r. Previous analyses have 
ignored limb-darkening and simply assumed a uniformly 
bright stellar disc with D — 1 for /i ^ = -^/l — {Ri,/rY 
and D = Q for fi < fii,. 

For a line with frequency integrated mass absorption 
coefficient k and of 'quality' q = vtht^/ (nac) ( Gayley||1995 1, 



we can write the perturbed direct component of the line ac- 
celeration in terms of an angle average of this core intensity 
(ORI, OP96) 

4:TVqKa 



-{^hD{r,fi) <56(r,^)). 



(1) 



where 6b is the perturbed escape probability. For an opti- 
cally thick line with Sobolev optical depth t^^ = KpL^ ^ 1, 
where = Vth/{dVn/dn) is the Sobolev length in direction 
n, we find for perturbations on length scales below that 
5b/5v oc fi/r^. Upon averaging, this leads to a strong insta- 
bility with growth rate Sgmd/Sv ~ v/L\. Since this is a fac- 
tor v/vth larger than the wind expansion rate dv/dr, small 
initial velocity perturbations are strongly amplified within 
this linear theory, by « Uoo/^^th ~ 100 e-folds (ORI). This 
implies such small-scale perturbations will quickly reach 
non-linear amplitudes within pure-absorption, non-Sobolev 



wind models, as first demonstrated by Owocki et al. ( 1988 1 



However, this strong de-shadowing instability can be 



counteracted by a 'line-drag' effect ( Lucy 1984 1 associated 
with the force of the diffuse, scattered radiation field. Ne- 
glecting perturbations in the source function S, the per- 
turbed diffuse force term is (ORII, OP96) 



c 



(2) 



where the minus sign signals the tendency of the line-drag 
to counteract the LDI. The similarity between eqs. [T] and [2] 
now allows us to write a very simple expression for the net 
relative reduction of the pure-absorption instability growth 
rate by the damping effect of this line-drag, 



ilnct{r) = ; = 1 



S{r){^^Sb{r,^I)} 



(3) 



A key issue for evaluating eqn. |3] lies in the computation of 
the assumed smooth source function S for the unperturbed 
background ffow. For a pure scattering line in a supersonic 
steady-state wind, this can be well approximated by the lo- 
cal Sobolev form 

;/*D(r,/i)6sob(r,/i)) 



(4) 



where the Sobolev escape probability 6sob('', m) = (1 ~ 
e~^f )/r^. For an optically thin line and a uniformly bright 
disc, the source function simply follows the dilution factor, 
S{r)/If: — (1— ^*)/2. (Eqn. A8 gives a corresponding expres- 
sion including limb-darkening.) For a thick line, this diluted 
form has a further correction for the angle-averaged escape 
probability (see OP96). 



2.1 Full angle integration 

For an optically thick line, applying eqn. [4] in eqn. [3] yields 
the asymptotic value finct — >■ 0.8 as r — >■ oo (ORII). Such 
a modest 20 % reduction of the pure-absorption instability 
growth rate implies the wind still is extremely unstable in 
its outer parts. 

But closer to the star the damping is stronger; in- 
deed, for a uniformly bright stellar disc with S{Ri,) = /*/2, 
the line-drag exactly cancels the de-shadowing instability 
at the stellar surface, i.e. Qnct — > when r — >■ 7i*. This 
holds for both optically thick and thin scattering source 
functions, and leads to marginal stabilization of the wind 
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base and to a later onset of wind structure than typically 
indicated by observations. However, noting that this ex- 
act cancellation occurs only because S{Ri,) — /*/2, let us 
now examine the efltect of photospheric limb-darkening on 
this net growth rate. For a simple Eddington (linear) limb- 
darkening law (eqn. A4|, analytic evaluation of eqn.js] gives 
f2net(i?*) ~ 0.06 and Q,nct{R*) ~ 0.22 for, respectively, an 
optically thick and thin source function. 

This means that in such limb-darkened models the line- 
drag no longer exactly cancels the LDI at the stellar surface 
and thus that also the wind base is unstable. 



2.2 Two-stream approximation 

In a numerical wind model aiming to simulate the non- 
linear evolution of the LDI, it is very computationally ex- 
pensive to carry out elaborate angle integrations at each 
radial grid-pint at each time-step. But testing has shown 
that a simple one-ray quadrature using a ray parameter 
y = {p/Ri,Y = 0.5, with impact parameter p = r^Jl — /i^, 
actually approximates rather well the full-angle integrated 
line force (see Appendix A), as long as separate accounts 
are taken for inward and outward rays in the diffuse force 
computations and a sphericity correction factor is applied 
(OP96). For this two-stream approximation, we find 



Table 1. Summary of stellar and wind parameters 



r2not(r) = 1 



(5) 



where fiy = yl — y(R*/r)^ is the local direction cosine of 
the ray at radius r. Assuming an optically thin source func- 
tion, this again results in a zero growth rate at the stellar 
surface for a uniform disc, for which eqn. [5| actually sim- 
plifies to Qnct{r) = A'*/(l + M*) (see also [Owocki fc Puis 



19991. By contrast, for a limb-darkened disc (again with an 
optically thin source function), we find $!(/?*) = 0.15, which 
is intermediate between the thin and thick results for full 
angle integration. 

Fig. [l] plots 0,nct{r) over r/Ri, = 1 — 1.5 for the cases 
discussed above, assu ming a smooth veloci ty law v = {1 — 
R^/rY with (3 = 0.8 ( |Pauldrach et al.|l986| ). Note that over 
this inner wind range, the two-stream approximation can 
again be viewed as an intermediate case between the fully 
angle-integrated cases with optically thick and thin source 
functions. 

Most significantly, since the pure-absorption growth 
rate Sgdir is of order 100 times the wind expansion rate, 
even the base- value of ~ 20 % of this rate found in this pa- 
per for limb-darkened models implies an absolute growth 
that is still substantially faster, by a factor ~ 20, than the 
wind expansion. This suggests that non-linear structure can 
now develop close to the wind base, as we next demonstrate. 



3 NUMERICAL SIMULATIONS OF THE 
TIME-DEPENDENT WIND 

To follow the non-linear evolution of such instability- 
generated wind structure, we must numerically integrate the 
radiation-hydrodynamical conservation equations of mass, 
momentum, and energy. Table [l] summarizes the stellar and 
wind parameters adopted in the models, which correspond 
to a typical early 0-supergiant in the Galaxy. The basic 



Name 


Parameter 


Value 


"^tf^llnT" 1 11 TTi 1 T^nci^'^r 

kj IjkjliCLl i lillilllUQl L y 


Li, 


S V 1 T,r^ 


Stellar mass 


Mi, 


50 Mq 


Stellar radius 


R* 


20 Rq 


Wind floor 






temperature 


7w 


40 000 K 


Initial steady-state 






- terminal speed 


Vao 


2000 km/s 


- mass-loss rate 


M 


2.1 X lO-^Mg/yr 


CAK exponent 


a 


0.65 


Line-strength 






- normalization 


Q 


2000 


- cut-off 


Qmax 


O.OIQ 


Ratio of ion tiiermal 






speed to sound speed 




0.28 


Electron scattering 






opacity 




0.34cm2/g 



simulation method used here has been extensively described 
in [Owocki eFaL] ( [19881 ), OP96, [Ruiiacres fc Owocki] ( [20021 ) ; 
this section recapitulates key assumptions and describes new 
features. 



3.1 Radiative cooling and radiative driving 

We solve the spherically symmetric conservation equations 
using the numerical hydrodynamics code VH-1 (developed 
by J. Blondin et al.). 

Radiative cooling of shock-heated gas is accounted for 
in the energy equation following Runacrcs & Owocki (2002j). 
This method mimics the effect of photoionization heating 
from the star's intense UV radiation field by simply de- 
manding that the gas never cools below a certain floor tem- 
perature, set here to 40 000 K (roughly the star's effective 
temperature) . 

A central challenge in these simulations is to compute 
the line-driving in a highly structured, time-dependent wind 
with a non-monotonic velocity. This requires a non-local in- 
tegration of the line-transport within each time-step of the 
simulation, in order to capture the instability near and be- 
low the Sobolev length. To allow for this objective, we adopt 



the smooth source function (SSF, Owocki 1991a I method 
described extensively in OP96. 

As in CAK, we assume a line-number distribution that 
is a power-law in line-strength q, but now (as in OP96) ex- 
ponentially truncated at a maximum strength Qmax, 



dN 
dq 



r{a)Q \Q 



-Q/Qn 



(6) 



with r(Qf) the complete gamma function. Here a is the CAK 
power-law index, which can be physically interpreted as the 
ratio of the line force due to optically thick lines to the total 
line force, and Q is a line-strength normalization constant, 
which can be interpreted as the ratio of the total line force 
to the electron scattering force in the case that all lines were 
optically thirj^ For typical O-supergiant conditions at solar 

^ Note that we have recast the line force using the Q notation of 
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metallicity, Qn 



Q ^ 2000 (Gayley 1995 Puis et al 



20001. In practice, keeping the nonlinear amplitude of the 



instabihty from exceeding the limitations of the numerical 



scheme requires a significantly smaller cut-off ( Owocki et al 



1988). First test-simulations (see also Runacres & Owocki 
2002 1 that increase Qmax to more physically realistic values 



indeed display stronger structure in wind regions where the 
instability is fully grown. But this higher Qmax also leads to 
spurious spikes in the velocity field, likely numerical arte- 
facts caused by the limiting grid resolution. The full impli- 
cations of the line-strength cut-off for wind structure is yet 
to be determined, and should be carefully examined in a 
detailed parameter study. However, we defer that to future 
work, since the test-simulations here also show that increas- 
ing the cut-off does not significantly affect the main focus of 
this paper, namely the onset of non-linear structure. 

All previous time-dependent wind simulations have as- 
sumed the stellar continuum radiation can be described by 
a uniformly bright stellar disc. But as discussed in ^ intro- 
ducing photospheric limb-darkening breaks the marginally 
stable properties of the wind base and so may lead to struc- 
ture formation also in the innermost wind. As also noted in 
^ for computational reasons we use a simple one-ray angle 
quadrature for both the direct and diffuse force components 
in all numerical models (see Appendix A) . Incorporating Ed- 
dington limb-darkening then typically increases the direct 
line force by a modest amount, ~ 3%, whereas the scat- 
tering source function (and thus the diffuse line force) at 
the stellar surface decreases by ~ 13%. This is the essential 
reason there is now a net instability at the wind base. 



3.2 Boundary conditions 

Each instability simulation evolves a smooth and relaxed 
Sobolev wind model. The lower boundary is set as in 
[Runacres fc Owocki] ( |2002[ ) , by fixing the density and tem- 
perature. Previous simulations have assumed the lowermost 
grid-point to be located at r = J?*, the photospheric radius 
where the continuum optical depth is unity, with a density 
po tuned to be approximately 5-10 times higher than at the 
sonic poinl[^ However, a simple ratio estimate of the steady- 
state sonic point density pa = M / {iivar^) to the typical pho- 
tospheric density p« ~ l/{HKe), with scale height H = /g, 
yields for the stellar and wind parameters in table [l] a much 
lower value, Pi^/p* < 0.01. This order-of-magnitude estimate 
is confirmed by more detailed, line-blanketed NLTE calcu- 
lations of steady-state, unified (photosphere-|-wind) model 
atmospheres of typical O supergiants like ( Pup (using the 
FASTWIND code. Puis et aLl|2005 |. In such models, we find 
that the sonic point is located ~ 1-2 % above the stellar sur- 
face, defined by ii* — r(TRoss = 2/3), and that the corre- 
sponding density contrast is p^j p-,, ~ 0.01. 



Gayley] l (l995l l rather than the kq notation of OP96. Q has the 
advantage of being a dimensionless measure of line-strength that 
is independent of the thermal speed. The relation between the 
two parameter formulations is Q = Ko'"thr(c«)^^'^~"V('^'^o). 
^ Using a higher value leads to long-lived oscillations at the Lamb 
frequency; using lower values risk choking the sclf-rcgulated mass 
loss of the line-driven wind ( [Owocki et al.^^1988, ^Owocki fc Puis] 
19991. ' ' 



The time-dependent simulations reported here thus as- 
sume r ~ 1.01-R* at the sonic point, which tends to shift 
inward the stellar surface somewhat as compared to earlier 
models. And since the stability properties of the lower wind 
are so sensitive to the specific value of the source function, 
which close to the wind base decreases very rapidly with 
radii, accounting for this shift can actually be quite signifi- 
cant (particularly for lower-density winds) . This adds to the 
effect of limb-darkening in making the region near the sonic 
point have a net instability. 

Moreover, if one accounts also for the likelihood that the 
stellar photosphere is not perfectly steady, but rather itself 
has a level of variability, this can further seed the growth 
of unstable structure in the overlying wind. Some previous 
simulations have explicitly perturbed the lower boundary 
by either sound waves or turbulence (e.g. 



Feldmeier et al 



1997 1, and indeed such models tend to show structure some- 



what closer to the wind base as compared to models with 



self-excited structure (compare Runacres & Owocki 2002 
with [Puis et al.||l993| and [Sundgvist et al.||2011| ). To study 
the influence of such perturbations in the presence of limb- 
darkening, the simulations here introduce a lower boundary 
sound wave of amplitude 5p/ po — 0.1 and period 4 000 s. 



4 NUMERICAL SIMULATION RESULTS 
4.1 Basic structure properties 

Let us first review the overall properties of numerical sim- 
ulations that follow the non-linear evolution of the compe- 



tition between the LDI and the line-drag (Owocki & Puis 
119991 IRunacres & Owocki||2002| IDessart & Owocki||2003|). 



The simulation displayed in Fig. [2| illustrates the formation 
of high-speed rarefactions that steepen into strong reverse 
shocks, which compress most of the wind material into dense 
and spatially narrow 'clumps' (or really shells in these 1-D 
simulations). This characteristic structure can be alterna- 
tively illustrated by plotting density and velocity against 
a Lagrangian mass coordinate tracking individual fluid ele- 
ments (as deflned by [Owocki fc Puls|1999[ their eqn. 14); the 
rightmost panel of Fig. |2] clearly shows how the high-speed 
flow consists of very rarefled gas and how most of the wind 
mass is concentrated into dense clumps. 

This model is computed assuming a uniformly bright 
stellar disc and without explicitly perturbing the lower 
boundary. The structure that evolves from the smooth wind 
at t = in the left panel of Fig.[2|is "self-excited", and arises 
from back-scattering of radiation from outer wind structure 
seeding small variations closer to the wind base, which are 



then subsequently amplified by the instability (Owocki & 
[Puls|1 999). 

Note how the line-drag effect discussed in ^ greatly 
suppresses instability growth close to the wind base; indeed, 
structure starts to develop only at r ~ 1.5i?* in this simula- 
tion, at odds with observations which typically indicate the 
presence of dense clumps much closer to the stellar surface 
(see The next section examines to what extent includ- 
ing limb-darkening and photospheric perturbations into the 
simulations can induce wind structure also at such low radii. 
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Figure 2. The left panel shows density and velocity contour plots of the time-evolution of a wind model with self-excited structure 
assuming a uniformly bright stellar surface. The middle and right panels show density and velocity of a single snapshot from the same 
model, taken ^ ISOksec after initiation, as function of radius (middle) and a Lagrangian mass coordinate (right). The dashed red line is 
the smooth start model. 




Figure 3. Simulated clumping factors /cK^). The left panel focuses on the inner wind regions r/Ri, = 1 — 2, whereas the right panel 
shows the wind all the way out to r/Ri, = 10 on a logarithmic abscissa. Blue dashed-dotted lines compare the simulations to the clumping 
law inferred for Pup by|Najarro et aL]||2011[|, see text. 'LD' indicates whether limb-darkening is accounted for in the models. 



V [km/s] log p [g/cm''] v [km/s] log p [g/cm'] 




Figure 4. Inner wind time evolutions of a simulation without limb-darkening and photospheric perturbations (left) and one including 
both effects (right). 
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4.2 Clumping in the inner wind 

Let us characterize the wind density structure in terms of a 
clumping factor 



(7) 



where the angle brackets here denote time averaging. Fig. [3] 
plots the radial variation of such clumping factors for the 
four simulation cases considered, i.e. for uniform-disc and 
limb-darkened models with and without explicit base pertur- 
bations, computed by averaging over each simulation time- 
step between t — 100 — SOOksec, where the simulations have 
become insensitive to the initial conditions. 

The right panel shows /d over the full simulation range 
r/7?* = 1 — 10, whereas the left panel focuses on the crucial 
inner wind regions where limb-darkening and base pertur- 
bations alter the onset of wind structure. The blue dashed 
curve in both panels show, for comparison, the empirically 
inferred clumping factoij^ for Pup based on the compre- 
hensive multi-diagnostic study by Najarro et al. (20111. 

The left panel illustrates that both perturbations and 
limb-darkening shift inward the onset of clumping as com- 
pared to the standard model. Much as anticipated in the 
linear stability analysis, limb-darkening increases the net 
growth rate, while the perturbation provides a seed which 
the instability can amplify to form the clumped structure. 
The right panel shows that this stronger structure persists 
to several stellar radii from the wind base. But note that 
our simulation volume here only extends to « lORi,; the 



analysis of Runacres & Owocki ( 2002 1 suggests that such 



strong structure will eventually dissipate and tend to settle 
at /cl « 4 in the outermost, radio emitting wind. 

We emphasize here that this initial exploration study 
only attempts to outline the general importance of limb- 
darkening and perturbations for the onset of LDI generated 
wind structure. Thus we do not aim for a perfect match 
to the empirical ^ Pup clumping law throughout the wind; 
quantitative details in the predicted clumping factors will 
depend on the exact treatment of limb-darkening and the 
source function, the nature and strength of the applied 
base perturbations, and the level of lateral fragmentation 
of clumps in multi-dimensional LDI simulations ( Dessart fc| 
|Owocki|2003[[2005a| ), as further discussed in !|5] Nonetheless, 
the new overall agreement between simulations and obser- 
vations in the inner wind is encouraging, and suggests that 
the standard line-deshadowing instability (possibly seeded 
by some modest photospheric perturbations) is indeed fully 
capable of generating significant structure also near the wind 
base. 



^ In such empirical work, /d is typically defined using volume- 
averaging rather than time-averaging, but in the stochastic 
medium here the two should be effectively interchangeable. If 
one n eglects the small inter-c lump density (but s ee |Zsarg6 et al.| 



2008 



Sundqvist et al. 



2010 



Surlan et al. 



20121, this clumping 



4.3 Time-evolution with and without 

limb-darkening and base perturbations 

Fig. |4] compares the time evolution of velocity and density 
in the inner wind for the model including limb-darkening 
and base perturbations (right) against the model ignoring 
both effects (left). Following the brief adjustment to initial 
conditions, the unperturbed model settles to a state with 
onset of intrinsic structure and variability at r ~ 1.5i?*. In 
contrast, the perturbed, limb-darkened model shows a much 
earlier onset that varies from a maximum r « 1.15J?* all the 
way down to the photospheric wind base. The characteristic 
time-scale of this variation, ~50ksec, is much longer than 
the perturbation period of ~4ksec. 

Our preliminary analysis suggests that this long-term 
variation is likely related to the intrinsic variations found 



in the "pure-absorption" models of Poe et al. (19901. As 



discussed there, these variations stem from a degeneracy of 
the CAK-like steady-state solutions in the regions near the 



wind sonic point. As shown in Fig. 7b of Poe et al. (19901, 



factor is equal to the inverse of the clump volume filling factor, 

/cl = l//vol- 



this leads to a wind-speed oscillation with period ~ 50 ksec. 
For slightly different conditions, the oscillating flow can ac- 
tually even stagnate and re-accrete onto the star (see Fig. 1 
in |Owocki|199lbl ). 

Remarkably, in the standard SSF model the introduc- 
tion of the diffuse force component eliminates this solution 
degeneracy, and so suppresses the oscillatory behaviour. But 
the additional introduction here of limb-darkening now en- 
hances the direct pure-absorption component of the line- 
force, and reduces the stabilizing diffuse component. This ev- 
idently allows the degeneracy to reappear, and so again leads 
to slow variations in the wind solution. The full consequences 
of this effect for wind structure and variability (for exam- 
ple, as a possible origin of discrete absorption components, 
Howarth fc Prinja|1989 1 must await multi-dimensional LDI 
wind models ( Dessart fc Owocki|2003 2005a I. Similarly, im- 
plications for wind initiation must await a future, more com- 
plete analysis of the dynamical role of diffuse radiation near 
the sub-sonic wind base (see, e.g., lOwocki fc Puls^^l999| . 
Note though, that in the current models the mass flux initi- 
ated at the wind base is still quite close to the initial steady 
state derived from standard CAK Sobolev theory. 



5 DISCUSSION, CONCLUSIONS, AND 
FUTURE WORK 

The central result of this paper is that photospheric limb- 
darkening in instability wind models leads to structure 
growth closer to the wind base than in previous models 
assuming a uniformly bright stellar disc. We demonstrate 
this both analytically in a linear perturbation analysis, and 
numerically using self-consistent, time-dependent radiation- 
hydrodynamical wind simulations. Particularly when com- 
bined with perturbations of the lower boundary, such limb- 
darkened models reproduce well the early onset of clumping 
typically inferred from observations of O-star winds. 

The analysis here uses a simple sound wave to exam- 
ine the general effects of lower boundary perturbations for 
the onset of LDI-generated wind structure. Physically, pho- 
tospheric perturbations could, for example, originate from 
non-radial pulsations, or even from the thin sub-surface con- 
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vective layer associated with the iron opacity peak ( Cantiello 
et al.|2009| |. 

Moreover, as in previous models, the mass flux at the 
wind base in these time-dependent simulations follows quite 
closely that of the initial steady state derived from standard 
CAK theory. This suggests the early onset of wind struc- 
ture found here does not directly alter the average mass-loss 
rate of the star. But note in this respect that the resulting 
clumped structure will modify e.g. the wind ionization bal- 
ance, and so indirectly affect the line-force. Such feedback ef- 
fects are not included in the simulations presented here, but 
may impact theoretical mass-loss rates derived from steady- 
state models fMuijres eFaLllWri ). 

For simplicity the simulations here assume spherical 
symmetry; first 2-D models by Dessart & Owocki (2003 



2005a I suggest somewhat lower clumping factors than in 



comparable 1-D models, by a factor of approximately two. 
Future work should carry out such multi-D LDI simulations 
including also the effects of limb-darkening and base pertur- 
bations. 

Finally, an unanticipated result of this study is the re- 
appearance of the solution degeneracy and associated slow 
global wind variations found in the pure-absorption mod- 
els of Poe et al. (19901. This arises from the change in the 



relative strength of the diffuse vs. direct force components 
in the transonic wind region of limb-darkened simulations. 
Determining the broad implications of these effects for wind 
structure and variability will require a careful analysis of 
the diffuse radiation in this wind initiation region around 
the sonic point. Thus, although standard CAK wind mod- 
els focus solely on the direct radiation from the star, this 
again (see also, e.g., Owocki & Puis 19991 emphasizes the 



subtle role of diffuse scattered radiation in the dynamics of 
line-driven winds. 
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APPENDIX A: THE EFFECT OF STELLAR 
LIMB-DARKENING ON THE RADIATION 
LINE FORCE 

Al The direct force term 

The direct line force term from an ensemble of lines at a 
distance r from the stellar center in spherical symmetry may 
be written as 



5dir(r) = g. 



- {fiD{p,r)b{^,r)) 



(Al) 



which is just OP96's eqn. 58 recast in Gayley s (1995) Q 
notation. Here QeQ = ffthin 
optically thin and 



is the line force if all lines were 



b{fi,r) 



(t>[x - fiv/vth 1 



[tQ{x,fJ',r, Zb) +Q/Qn 



-dx 



(A2) 



the ensemble-summed escape probability (OP96, eqn. 57), 
where is the profile-function (assumed here to be a Gaus- 
sian) and tQ{x,^,r,Zh) the frequency-dependent optical 
depth for a line of strength q = Q. 

Assuming a uniformly bright stellar disc, and casting 
the angle integration in terms of a ra y parameter y 
{p/Ri,)^ (see |2|, we may write eqn. 



Al 



ffdir(j-) 



Sthi 



b{fiy,r)dy, 



(A3) 



with /iy(r) = ^yi — y{Ri,/r)^ the direction cosine of the ray. 

Let us next consider a linear Eddington limb-darkening 
law (see, e.g., Mihalas|1978 eqn. 3.37), 



D(M,r) 



13,13 
2 + 4^ ^2 + 4 



1 3 n 

2 + 4V^' 



(A4) 



where ji = cos 9' is the direction cosine between a star- 
centered radius vector and the ray at the stellar surface, 
and where /^t and /i are related to, respectively, the dilution 
factor and impact parameter (as defined in ipk. 

Applying this limb-darkening law in eqn. [A1| yields 



Id / N 



Pthin X 



6(/iy,r-)dy + +- / ^l~yb{^iy,r)dy\ . (A5) 



In a time-dependent wind model, these angle integrals must 
be evaluated numerically. However, to carry out elaborate 
angle integrations at each radial grid-pint at each time- 
step is very computationally expensive. Fig. |A1| plots line 
forces computed at the first evolution time-step of a previ- 
ously relaxed Sobolev model, using i) full numerical angle- 
integrations of eqn. A3 and A5 and ii) a simple one-ray 
quadrature with y = 0.5. The figure shows that for both 
limb-darkened and uniform-disc models this one-ray quadra- 
ture gives quite accurate results, within 10-20 % of the full 
angle integration in the sub-sonic region and much better 
in the part of the wind where structure develops. As in pre- 
vious uniform-disc LDI simulations (e.g. 



Puis 1999 Runacres & Owocki 20021, all time-dependent 



OP96; Owocki & 



wind models presented in the main text of this paper use 
this one-ray formulation. 
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Figure Al. Upper panels show ratios of direct line- forces com- 
puted with and without limb-darkening (LD), for full angle- 
integrations (FI) and a one-ray quadrature with y = 0.5. Lower 
panels show ratios of forces computed with full angle-integrations 
and the same one-ray quadrature, for uniform-disc and limb- 
darkened models. The right panels focuses on the innermost wind, 
with the blue dotted lines marking the position of the steady-state 
sonic point. 



Note further from Fig. |Al| that nowhere in the wind is 
the effect of limb-darkening on the direct line force bigger 
than ~ 10 %. 



A2 The diffuse force term 

The diffuse line force in the smooth source function (SSF) 
approximation is (OP96, eqn 62) 



ffdiff(r) = — 2gthin 



Sir) 



[fib(^j,,r)). 



(A6) 



Using again the ray parameter y, we may approximate 
eqn. |A6| by 



5diff(?-) 



Sir) 



Stir) 2(1- 
1 



(6_(/iy,r) - b+iiiy,r))dy, 



(A7) 



where we have now taken separate accounts for outward 
streaming ("+", fi > 0) and inward streaming ("— ", pi < 0) 
photons, and applied a simple (r/i?*)^ correction factor to 
account for the fact that the angle integ ral here really should 
extend to y = ( |Owocki||l991a OP96). Expressions 

for the corresponding optical depths can be found in OP96. 

5't(r)//* = (1 — /i*)/2 is here the optically thin, pure 
scattering source function in the case of a uniformly bright 
stellar disc (eqn.|4|. Applying the Eddington limb-darkening 
law eqn. |A4| in eqn. [4] this optically thin source function now 
takes the form 



Sl^ir) 



, ln[/i,/(l + Vl - /i?)! \ 
3t4 /— ^ . (A8) 



which takes the value 7/16 (i.e. < 1/2) at the stellar surface, 
where /i* — > 0. By simply applying this source function in 
eqn. |A7| we can then approximate the effect of stellar limb- 
darkening on the diffuse line-force component. 



Figure A2. Finite-disc correction factors for a uniform (black) 
and limb-darkened (red) disc, assuming a smooth /3 = 0.8 velocity 
law and a = 0.65. 



A3 Sobolev approximation 

For completeness, we add corresponding limb-darkening 
modifications also for the Sobolev line force. Now eqn. |A2| 
can be solved analytically to yield for the CAK direct line 
force (which is also the total line force since in the Sobolev 
approximation the diffuse force vanishes everywhere) 



Sob 
5dir 



5tliin 



(1 - a)i 



(A9) 



with Sobolev optical depths r = Qt = Qk^pc/ idv /dr) and 



fir) 



gSob 
„rad 
ySob 



t_t, and finite-disc correction factor 
2 



(l-M?)(l + a)- 
1 



(AlO) 



DitJ,,r)^il + ati^)°'dfi, (All) 
; d(ln ?;)/d(ln r) — 1 is the so-called wind anisotropy 

1, 



(A12) 



where a ; 
factor. 

For the case of a uniformly bright disc with D 
eqn. |A11| becomes 

fM ^ il + a)'+"-il + a^^l)'+- 
^ il + a)a il + a)'-il - fii) 

As r — >■ oo, |J,^, 1 and a — >■ —1, and thus / — ^ 1, as 
expected since at large distances the star closely resembles 
a point source. But closer to the star, where cr > 1, the 
finite-disc correction reduces the line-force somewhat, with 
the limiting value /(r — > 7?*) — > 1/(1 -I- a). 

For a limb-darkened disc, eqn. |A11| cannot generally 
be solved analytically. However, again taking the near-star 



limit for the simple Eddington-case, one obtains (Cranmer 
|fc Owocki|1995t 



fiR*) = 



1 3 

2 2 



1 + Q 
3 + 2a 



(A13) 



The square brackets are here the limb-darkening correction 
factor, equal to unity for optically thin lines with q = 0, but 
always somewhat greater than unity (up to 1.2 for a = 1) 
for a realistic line-ensemble with optically thick line-fraction 
a > 0, simply due to the more centrally concentrated stellar 
light. 
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Fig. |A2| compares the two radial finite-disc functions for 
Q = 0.65 and a canonical v = Uoo(l — Ri,/rY velocity law, 
using P — 0.8. This again shows that the general effect of 
limb-darkening is modest, with the minimum ratio //./''^ ~ 
0.93 occurring at the stellar surface. Such quite small effects 
are in accordance with the results found in the previous 
section for non-Sobolev line forces. 

Indeed, from numerical hydrodynamics simulations that 
relax a Sobolev line-force model, including full angle- 
integrations of the limb-darkened finite-disc correction 
throughout the wind, we find only a modest « 10% higher 
mass-loss rate than in comparable models that assume a 
uniformly bright stellar disc. 

As demonstrated below, this small effect can be read- 
ily understood analytically, even though the complex spatial 
and velocity dependence of the finite-disc correction factor 
generally prohibits full analytic solutions for the steady- 
state equation of motion. However, for the uniformly bright 
stellar-disc case, full numerical solutions (e.g.. Friend & Ab- 
bott||1986 Pauldrach et al.|1986 l show a reduced mass-loss 
by a factor « f{R^,y'°' = 1/(1 -I- a)^/" as compared to the 
original CAK point star rate, 



M 



L 



l/o, 



l-Fc 



(A14) 



where a line-strength cut-off Tmax 2> 1 has been assumed. 

We may understand this effect of the finite extent of 
the stellar disc on the mass- loss rate by noting that, since 
/(r) increases outward, the stellar surface now represents a 
"critical point" from which it is hardest to accelerate the 
wind, thereby also setting the maximal allowed mass-loss 
rate. Thus the simple substitution / — >■ shows that in- 
cluding photospheric limb-darkening increases the finite-disc 
corrected mass-loss rate by « 10% for reasonable values 
Q ~ 2/3 (see also Cranmer & Owocki 1995 Cure et al 



2012), as also found above for the numerical model. 
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